Semivariogram Clouds and Plots

Outline

The dataset we are using is from the data paper https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/14-1345.1 by Stallings et al. We will use salinity data collected at 80 sites from the Big Bend down to the Springs Coast of Florida. In the code chunks below we will cover

• Reading and Writing Spatial Data
• Summarizing Distances
• Clouds and Semivariograms
• Randomization of Semivariograms

Libraries and references

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.3     v purrr   0.3.4
## v tibble  3.1.1     v dplyr   1.0.5
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(sf)

## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

# edit May 28, 2021 - install mapview 2.9.9 from github and then knitting to html will work
remotes::install_github("r-spatial/mapview")

## Downloading GitHub repo r-spatial/mapview@HEAD

## processx    (3.5.1  -> 3.5.2 ) [CRAN]
## e1071       (1.7-6  -> 1.7-7 ) [CRAN]
## rlang       (0.4.10 -> 0.4.11) [CRAN]
## colorspace  (2.0-0  -> 2.0-1 ) [CRAN]
## systemfonts (1.0.1  -> 1.0.2 ) [CRAN]
## vctrs       (0.3.7  -> 0.3.8 ) [CRAN]
## pillar      (1.6.0  -> 1.6.1 ) [CRAN]
## fansi       (0.4.2  -> 0.5.0 ) [CRAN]
## ellipsis    (0.3.1  -> 0.3.2 ) [CRAN]
## tibble      (3.1.1  -> 3.1.2 ) [CRAN]
## xfun        (0.22   -> 0.23  ) [CRAN]
## viridis     (0.6.0  -> 0.6.1 ) [CRAN]
## raster      (3.4-5  -> 3.4-10) [CRAN]
## sfheaders   (NA     -> 0.4.0 ) [CRAN]
## rapidjsonr  (NA     -> 1.2.0 ) [CRAN]
## jsonify     (NA     -> 1.2.1 ) [CRAN]
## geometries  (NA     -> 0.2.0 ) [CRAN]
## httpuv      (1.6.0  -> 1.6.1 ) [CRAN]
## geojsonsf   (NA     -> 2.0.1 ) [CRAN]
## servr       (NA     -> 0.22  ) [CRAN]
## leafpop     (0.0.6  -> 0.1.0 ) [CRAN]
## leafem      (0.1.3  -> 0.1.6 ) [CRAN]

## Installing 22 packages: processx, e1071, rlang, colorspace, systemfonts, vctrs, pillar, fansi, ellipsis, tibble, xfun, viridis, raster, sfheaders, rapidjsonr, jsonify, geometries, httpuv, geojsonsf, servr, leafpop, leafem

## Warning: package 'tibble' is in use and will not be installed

## Installing packages into 'C:/R-4.0.5/library'
## (as 'lib' is unspecified)

## 
##   There are binary versions available but the source versions are later:
##         binary source needs_compilation
## e1071    1.7-6  1.7-7              TRUE
## fansi    0.4.2  0.5.0              TRUE
## leafpop  0.0.6  0.1.0             FALSE
## leafem   0.1.3  0.1.6             FALSE
## 
##   Binaries will be installed
## package 'processx' successfully unpacked and MD5 sums checked
## package 'e1071' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'e1071'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem
## copying C:\R-4.0.5\library\00LOCK\e1071\libs\x64\e1071.dll to C:
## \R-4.0.5\library\e1071\libs\x64\e1071.dll: Permission denied

## Warning: restored 'e1071'

## package 'rlang' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'rlang'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem
## copying C:\R-4.0.5\library\00LOCK\rlang\libs\x64\rlang.dll to C:
## \R-4.0.5\library\rlang\libs\x64\rlang.dll: Permission denied

## Warning: restored 'rlang'

## package 'colorspace' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'colorspace'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying
## C:\R-4.0.5\library\00LOCK\colorspace\libs\x64\colorspace.dll to C:
## \R-4.0.5\library\colorspace\libs\x64\colorspace.dll: Permission denied

## Warning: restored 'colorspace'

## package 'systemfonts' successfully unpacked and MD5 sums checked
## package 'vctrs' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'vctrs'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem
## copying C:\R-4.0.5\library\00LOCK\vctrs\libs\x64\vctrs.dll to C:
## \R-4.0.5\library\vctrs\libs\x64\vctrs.dll: Permission denied

## Warning: restored 'vctrs'

## package 'pillar' successfully unpacked and MD5 sums checked
## package 'fansi' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'fansi'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem
## copying C:\R-4.0.5\library\00LOCK\fansi\libs\x64\fansi.dll to C:
## \R-4.0.5\library\fansi\libs\x64\fansi.dll: Permission denied

## Warning: restored 'fansi'

## package 'ellipsis' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'ellipsis'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem
## copying C:\R-4.0.5\library\00LOCK\ellipsis\libs\x64\ellipsis.dll to C:
## \R-4.0.5\library\ellipsis\libs\x64\ellipsis.dll: Permission denied

## Warning: restored 'ellipsis'

## package 'xfun' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'xfun'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem
## copying C:\R-4.0.5\library\00LOCK\xfun\libs\x64\xfun.dll to C:
## \R-4.0.5\library\xfun\libs\x64\xfun.dll: Permission denied

## Warning: restored 'xfun'

## package 'viridis' successfully unpacked and MD5 sums checked
## package 'raster' successfully unpacked and MD5 sums checked
## package 'sfheaders' successfully unpacked and MD5 sums checked
## package 'rapidjsonr' successfully unpacked and MD5 sums checked
## package 'jsonify' successfully unpacked and MD5 sums checked
## package 'geometries' successfully unpacked and MD5 sums checked
## package 'httpuv' successfully unpacked and MD5 sums checked
## package 'geojsonsf' successfully unpacked and MD5 sums checked
## package 'servr' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\Meagan.Schrandt\AppData\Local\Temp\RtmpG2sRMj\downloaded_packages

## installing the source packages 'leafpop', 'leafem'

## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install Rtools 4.0 from https://cran.r-project.org/bin/windows/Rtools/.

##   
  
  
   checking for file 'C:\Users\Meagan.Schrandt\AppData\Local\Temp\RtmpG2sRMj\remotesb2471bc1c80\r-spatial-mapview-125842c/DESCRIPTION' ...
  
   checking for file 'C:\Users\Meagan.Schrandt\AppData\Local\Temp\RtmpG2sRMj\remotesb2471bc1c80\r-spatial-mapview-125842c/DESCRIPTION' ... 
  
v  checking for file 'C:\Users\Meagan.Schrandt\AppData\Local\Temp\RtmpG2sRMj\remotesb2471bc1c80\r-spatial-mapview-125842c/DESCRIPTION' (412ms)
## 
  
  
  
-  preparing 'mapview': (726ms)
##    checking DESCRIPTION meta-information ...
  
   checking DESCRIPTION meta-information ... 
  
v  checking DESCRIPTION meta-information
## 
  
  
  
-  checking for LF line-endings in source and make files and shell scripts (443ms)
## 
  
  
  
-  checking for empty or unneeded directories
## 
  
  
  
-  building 'mapview_2.9.9.tar.gz'
## 
  
   
## 

## Installing package into 'C:/R-4.0.5/library'
## (as 'lib' is unspecified)

library(mapview)
library(leaflet)
library(leafpop) # for popups in mapviews
library(gstat) # for semivariograms
library(lattice) # for random semivariogram plots
mapviewOptions(fgb = FALSE)
# setwd("E:/bigbend_grass")

# some spatial data analysis references
# https://zia207.github.io/geospatial-r-github.io/semivariogram-modeling.html
# https://pdixon.stat.iastate.edu/stat406/notes/part%203b-4.pdf

# https://community.esri.com/t5/arcgis-geostatistical-analyst/quot-spatial-statistical-data-analysis-for-gis-users-quot/td-p/394418

Reading in Data

This chunk takes a data frame with geographic coordinates and converts it to a simple features (sf) object that uses the project coordinate reference system of Albers Equal Area, which has the units of meters.

## 
## -- Column specification --------------------------------------------------------
## cols(
##   .default = col_double(),
##   Site = col_character(),
##   Date = col_character(),
##   Time = col_time(format = ""),
##   Notes = col_character()
## )
## i Use `spec()` for the full column specifications.

## cols(
##   xx = col_double(),
##   Site = col_character(),
##   Season = col_double(),
##   Date = col_character(),
##   Time = col_time(format = ""),
##   Depth = col_double(),
##   Towdist = col_double(),
##   Towarea = col_double(),
##   SGcont = col_double(),
##   perTha = col_double(),
##   perSyr = col_double(),
##   perOSG = col_double(),
##   perTotSG = col_double(),
##   lenTha = col_double(),
##   lenSyr = col_double(),
##   lenOSG = col_double(),
##   perSand = col_double(),
##   perAlgae = col_double(),
##   volAlgae = col_double(),
##   denAlgae = col_double(),
##   Secchi = col_double(),
##   Landist = col_double(),
##   Edgedist = col_double(),
##   `1*Lhab` = col_double(),
##   `2*Lhab` = col_double(),
##   Long = col_double(),
##   Lat = col_double(),
##   `Segment (BB vs SC)` = col_double(),
##   Salinity = col_double(),
##   Conduct = col_double(),
##   DO = col_double(),
##   DOsat = col_double(),
##   Temp = col_double(),
##   Notes = col_character()
## )

##  [1] "xx"                 "Site"               "Season"            
##  [4] "Date"               "Time"               "Depth"             
##  [7] "Towdist"            "Towarea"            "SGcont"            
## [10] "perTha"             "perSyr"             "perOSG"            
## [13] "perTotSG"           "lenTha"             "lenSyr"            
## [16] "lenOSG"             "perSand"            "perAlgae"          
## [19] "volAlgae"           "denAlgae"           "Secchi"            
## [22] "Landist"            "Edgedist"           "1*Lhab"            
## [25] "2*Lhab"             "Long"               "Lat"               
## [28] "Segment (BB vs SC)" "Salinity"           "Conduct"           
## [31] "DO"                 "DOsat"              "Temp"              
## [34] "Notes"

## spec_tbl_df[,34] [170 x 34] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ xx                : num [1:170] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Site              : chr [1:170] "BB09_001" "BB09_002" "BB09_003" "BB09_004" ...
##  $ Season            : num [1:170] 2009 2009 2009 2009 2009 ...
##  $ Date              : chr [1:170] "7/1/2009" "7/9/2009" "7/14/2009" "6/26/2009" ...
##  $ Time              : 'hms' num [1:170] 08:55:00 18:00:00 09:00:00 10:10:00 ...
##   ..- attr(*, "units")= chr "secs"
##  $ Depth             : num [1:170] 2.1 1.3 1.8 5.1 3.4 3.3 3.4 3.2 2.4 3.1 ...
##  $ Towdist           : num [1:170] 86 34 80 79 88 100 81 96 89 76 ...
##  $ Towarea           : num [1:170] 160.8 63.6 149.6 147.7 164.6 ...
##  $ SGcont            : num [1:170] 1 2 2 1 2 1 2 1 2 2 ...
##  $ perTha            : num [1:170] 0.1 0.85 0.25 0.1 0.25 0.6 0.9 0.9 0 0.25 ...
##  $ perSyr            : num [1:170] 0 0 0 0.1 0.75 0 0.05 0 0.2 0 ...
##  $ perOSG            : num [1:170] 0 0 0 0.05 0 0 0 0 0 0 ...
##  $ perTotSG          : num [1:170] 0.1 0.85 0.25 0.25 1 0.6 0.95 0.9 0.2 0.25 ...
##  $ lenTha            : num [1:170] 10 15 10 20 40 35 40 15 NA 25 ...
##  $ lenSyr            : num [1:170] NA NA NA 20 60 NA NA NA 35 NA ...
##  $ lenOSG            : num [1:170] NA NA NA 10 NA NA NA NA NA NA ...
##  $ perSand           : num [1:170] 0.9 0.15 0.2 0.65 0 0.4 0.05 0.1 0.8 0.15 ...
##  $ perAlgae          : num [1:170] 0 0 0.55 0.1 0 0 0 0 0 0.6 ...
##  $ volAlgae          : num [1:170] 0 15142 17034 1875 17034 ...
##  $ denAlgae          : num [1:170] 0 238.2 113.9 12.7 103.5 ...
##  $ Secchi            : num [1:170] 1.2 1.3 1.8 5.1 NA 3.3 3.4 3.2 1.8 3.1 ...
##  $ Landist           : num [1:170] 1.2 3.3 3.4 11.9 9.8 10.2 13.6 8.5 5.8 12.7 ...
##  $ Edgedist          : num [1:170] 0.3 13.9 12 0.3 3.2 5.8 2.2 8 0.5 3.5 ...
##  $ 1*Lhab            : num [1:170] 1 1 5 1 1 1 1 1 1 1 ...
##  $ 2*Lhab            : num [1:170] 2 5 1 5 2 5 2 5 5 5 ...
##  $ Long              : num [1:170] -83.2 -82.7 -82.7 -83.8 -82.8 ...
##  $ Lat               : num [1:170] 29.4 28.6 28.4 29.8 28.9 ...
##  $ Segment (BB vs SC): num [1:170] 0 1 1 0 1 0 1 0 0 1 ...
##  $ Salinity          : num [1:170] NA NA NA NA NA NA NA NA NA NA ...
##  $ Conduct           : num [1:170] NA NA NA NA NA NA NA NA NA NA ...
##  $ DO                : num [1:170] NA NA NA NA NA NA NA NA NA NA ...
##  $ DOsat             : num [1:170] NA NA NA NA NA NA NA NA NA NA ...
##  $ Temp              : num [1:170] NA 30.3 30.4 29 NA 28 30.6 29 NA 29 ...
##  $ Notes             : chr [1:170] "Thal 10-20cm in small patches" "Thal 10-20cm; towed seagrass \"island\" in macroalgae \"field\"" "Thal <= 10cm" "HAWR 10cm; Thal/Syr 10-30cm; temp from depth sounder" ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   xx = col_double(),
##   ..   Site = col_character(),
##   ..   Season = col_double(),
##   ..   Date = col_character(),
##   ..   Time = col_time(format = ""),
##   ..   Depth = col_double(),
##   ..   Towdist = col_double(),
##   ..   Towarea = col_double(),
##   ..   SGcont = col_double(),
##   ..   perTha = col_double(),
##   ..   perSyr = col_double(),
##   ..   perOSG = col_double(),
##   ..   perTotSG = col_double(),
##   ..   lenTha = col_double(),
##   ..   lenSyr = col_double(),
##   ..   lenOSG = col_double(),
##   ..   perSand = col_double(),
##   ..   perAlgae = col_double(),
##   ..   volAlgae = col_double(),
##   ..   denAlgae = col_double(),
##   ..   Secchi = col_double(),
##   ..   Landist = col_double(),
##   ..   Edgedist = col_double(),
##   ..   `1*Lhab` = col_double(),
##   ..   `2*Lhab` = col_double(),
##   ..   Long = col_double(),
##   ..   Lat = col_double(),
##   ..   `Segment (BB vs SC)` = col_double(),
##   ..   Salinity = col_double(),
##   ..   Conduct = col_double(),
##   ..   DO = col_double(),
##   ..   DOsat = col_double(),
##   ..   Temp = col_double(),
##   ..   Notes = col_character()
##   .. )

## [1] "sf"         "tbl_df"     "tbl"        "data.frame"

## Reading layer `bigbend_grass' from data source `C:\Users\Meagan.Schrandt\OneDrive - Florida Fish and Wildlife Conservation\Coding_General\R\FWRI_R_Club\Scripts_Mtg_Docs\2021-05-11_FWRI_RClub_spatial_variogram\bigbend_grass.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1902 features and 8 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 381734.7 ymin: 461716.6 xmax: 538871.2 ymax: 678105.8
## Projected CRS: Albers Conical Equal Area (Florida Geographic Data Library)

## [1] FALSE

## [1] TRUE

Writing sf objects to shapefiles and geopackage geodatabase

A sf object can be written out as an ESRI shapefile. If you have several sf objects, such as points for sites and polygons of seagrass areas, you can keep them together by writing them out as a geopackage geodatabase.

# Writing a single shapefile to the working directory
st_write(hab_sf1,
dsn = "hab1.shp",
driver = "ESRI Shapefile",
 append = FALSE)

# This will put the sf layer into the geodatabase
st_write(hab_sf1, dsn = file.path(getwd(), "bbsg_v1_2021_03_18.gpkg"), layer = "habitat", driver = "GPKG", quiet = FALSE)

# We can check on what layers are in the geodatabase.
st_layers("bbsg_v1_2021_03_18.gpkg") 
 
# Now we can put in the polygons.
 st_write(bigbend_seagrass, dsn = file.path(getwd(), "bbsg_v1_2021_03_18.gpkg"), layer = "bigbend_seagrass", driver = "GPKG", quiet = FALSE, delete_layer=TRUE)
st_layers("bbsg_v1_2021_03_18.gpkg") 

# Note if we were to modify hab_sf1 in R, for example, then to replace the existing layer you have to specify both delete_layer = TRUE and append = TRUE arguments so that the old sf object is deleted and the new one is appended to replace it.

# st_write(hab_sf1, dsn = file.path(getwd(), "st_layers("bbsg_v1_2021_03_18.gpkg") "), layer = "hab_sf1", delete_layer = TRUE, append = TRUE, driver = "GPKG", quiet = FALSE)

Summary of Distances

Getting a summary of distances is important because we are going to be looking at how spatial variance, the semivariance of salinity, is a function of distance.

# use st_distance function to create pairwise distance matrix among all the sites
# extract half of that matrix and convert it to a vector then get the summary of the distances, min, Q1, median, mean, Q3, and max
hab_sf2 <- filter(hab_sf1, (!is.na(Salinity)))
distmat_hab <- st_distance(hab_sf2, hab_sf2,by_element = FALSE)

dim(distmat_hab) # get the dimensions of the distance matrix

## [1] 80 80

rdistmat_hab <- distmat_hab[1:80, ] # change to match # sites/rows
rdistmat1_hab <- as.vector(rdistmat_hab)

round(summary(rdistmat1_hab[rdistmat1_hab!=0]))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     807   27764   94038  104561  180636  243452

# Knowing max distance important for variogram function as that is used in the default cutoff. That cutoff is the spatial separation distance up to which a pair of points are included in semivariance estimates.  The default is the length of the diagonal of bounding box spanning the data, which is ~ approximately the max distance.  Our default cutoff will be 243,452/3 = 81,150 meters or ~ 82 km.

Clouds and Semivariograms

We will start with a semivariogram cloud to examine spatial variation in salinity. Specifically, just we had calculated all pairwise distances between sites, we will calculate all pairwise semivariances in salinity, and plot those semivariances as a function of their distances.

A semivariogram cloud presents a smear of points, but by binning those points into distance classes and by binning the semivariances of those points through the width argument we can make a semivariogram plot. With the semivariogram plot we are asking: do stations near each other have similar salinities? Is the semivariance, gamma, small at smaller distances and does it increase with larger distances?

ggplot(hab_sf2, aes(Salinity)) + geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# histogram of salinity

# log transformation not help so stick with original data
ggplot(hab_sf2, aes(log(Salinity))) + geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(hab_sf2$Salinity)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   14.00   23.50   28.50   27.22   29.95   34.70

# Semivariogram Cloud
z_sal_cloud =  variogram(Salinity ~ 1, hab_sf2, cloud = TRUE)
class(z_sal_cloud)

## [1] "variogramCloud" "data.frame"

plot(z_sal_cloud)

# quite of smear of points
z1 <- as.data.frame(z_sal_cloud)
View(z1) # brings up the semivariance cloud data so we can sort on gamma, the semivariance, to find our maximum gamma occurs at a pair of sites about 26 km apart. But, where are those sites and what are their salinities?

# https://r-spatial.github.io/gstat/reference/plot.variogramCloud.html

# the inner plot statement is saying plot the cloud as we did above, but with identiy = TRUE we now get a interactive crosshairs to select a point on graph.  Once we make that selection and hit finish, then the outer plot statement executes to map our sites and connect the pairs of sites with a red line
plot(plot(z_sal_cloud, identify = TRUE), hab_sf2)

## [1] "mouse-left identifies, mouse-right or Esc stops"

# now going back to cloud plot we see a pair of numbers identifying that point.  We can view z1 and enter 75 under left and 56 under right to filter to only that point.   
view(hab_sf2)
# 56 and 75 are the the row numbers in the data frame hab_sf2 so we can get get the salinity measurements and hand calculate gamma in z1.
# 56, site BB10_087, has a salinity of 16
# 75, site BB10_116, has a salinity of 34.7

# gamma = 1/2(z1 - z2)^2 = 1/2(34.7 - 16)^2 = 1/2(18.7^2) = 1/2(346.9)
# gamma = 174.845

# Semivariogram Plot
# accept the defaults
z_sal_sv1 = variogram(Salinity ~1, hab_sf2)
plot(z_sal_sv1, ylab = "Salinity Semivariance")

# note the smaller scale now on the y-axis
# how are these 15 points, and their semivariance and distance, related to the cloud we saw earlier?

plot(z_sal_sv1, plot.numbers = TRUE) #shows number of pairs in each binned point and we see that in the corresponding semivariogram object

View(z_sal_sv1)

# now specify cutoff and width arguments.  width specifies the distance size class, or bin or lag size, over which we will average distances and semivariances we had in the cloud.  The specifications below mean we look at those points from the cloud in 10 km increments out to 80 km so this semivariogram plot will have 8 points plotted.
z_sal_sv2 = variogram(Salinity ~1, hab_sf2, cutoff = 80000, width = 10000)
plot(z_sal_sv2, ylab = "Salinity Semivariance")

View(z_sal_sv2)
# now view cloud data frame and filter for distances (dist) 0-10000
View(z1)
# those are 254 cloud points that had both averages taken of distance and gamma to make that first point in the semivariogram plot

filter(z_sal_cloud, dist <= 10000) %>% 
  summarize(
    np = dplyr::n(),
    mean_dist = mean(dist),
    mean_gamma = mean(gamma)
  )

##    np mean_dist mean_gamma
## 1 254   6164.45    12.5002

# these results match the z_sal_sv2 results

Random semivariograms

I downloaded Chapter 8 code from Applied Spatial Data Analysis with R, https://asdar-book.org/, for generating random variograms. An initial empirical, or sample, variogram is calculated and plotted for salinity, and that object is saved. Then salinity values are randomly sampled and assigned to the coordinates so that a random semivariogram can be calculated. That is done 100 times, and the grey lines of the random semivariograms are plotted with the blue line of the empirical semivariogram. The semivariogram for salinity shows suggestive evidence of spatial autocorrelation out to ~ 10 km compared to randomized semivariograms. I have used the script below to explore spatial structure in National Aquatic Resource Surveys (NARS) data.

print(xyplot(gamma ~ dist, z_sal_sv1, pch = 3, type = 'b', lwd = 2, col = 'darkblue',
             panel = function(x, y, ...) {
               for (i in 1:100) {
                 hab_sf2$random = sample(hab_sf2$Salinity)
                 v = variogram(random ~ 1, hab_sf2)
                 llines(v$dist, v$gamma, col = 'grey')
               }
               panel.xyplot(x, y, ...)
             },
             ylim = c(0, 35), xlab = 'distance', ylab = 'Salinity semivariance'
))

# Semivariograms are sensitive to outliers.  See https://pdixon.stat.iastate.edu/stat406/notes/part 3b-4.pdf
# so use Cressie-Hawkins to create robust variogram

z_sal_sv3 = variogram(Salinity ~1, hab_sf2, cressie = TRUE)
plot(z_sal_sv3, cressie =TRUE, ylab = "Salinity Robust Semivariance")

# change ylim below in random semivariogram based on this plot

print(xyplot(gamma ~ dist, z_sal_sv3, pch = 3, type = 'b', lwd = 2, col = 'darkblue',
             panel = function(x, y, ...) {
               for (i in 1:100) {
                 hab_sf2$random = sample(hab_sf2$Salinity)
                 v = variogram(random ~ 1, hab_sf2, cressie = TRUE)
                 llines(v$dist, v$gamma, col = 'grey')
               }
               panel.xyplot(x, y, ...)
             },
             ylim = c(0, 35), xlab = 'distance', ylab = 'Salinity robust semivariance'
))

The marked difference in the first two bins of classical semivariance versus Cressie-Hawkins semivariance reflects what we saw in the semivariance cloud. We have pairs of sites at close distances, < 25 km, having large semivariance values. In other words, low salinity sites close to high salinity sites.

Disclaimer

The United States Environmental Protection Agency (EPA) project code is provided on an “as is” basis and the user assumes responsibility for its use. EPA has relinquished control of the information and no longer has responsibility to protect the integrity , confidentiality, or availability of the information. Any reference to specific commercial products, processes, or services by service mark, trademark, manufacturer, or otherwise, does not constitute or imply their endorsement, recommendation or favoring by EPA. The EPA seal and logo shall not be used in any manner to imply endorsement of any commercial product or activity by EPA or the United States Government.